package Math::Polynomial::RealRoots;
use strict;
use Carp;

use Exporter 'import';
our @EXPORT_OK = qw(real_roots);
use constant {
	DBG	=> 0,
	EPSILON	=> 1e-16,
};

sub diff(@) {
	my @d;
	for (my $i = $#_; $i > 0; $i--) {
		$d[$i-1] = $_[$i] * $i
	}
	@d
}

sub real_roots(@) {
	my @coef = @_;
	my @roots;

	pop @coef while @coef && !$coef[-1];
	return poly_solve(@coef);
}

# util. evaluate poly at value, not used by solver proper
sub peval {
	my ($p, $x) = @_;
	my $s = $p->[-1];
	for (my $i = $#$p-1; $i >= 0; $i--) {
		$s = $s * $x + $p->[$i]
	}
	$s
}

# calculates a range [-max, max] such that there's no
# root outside of it
sub upper_bound {
	my ($nn, $vv, $max) = ($#_, $_[-1], -1e60);
	for (my $i = $nn - 1; $i >= 0; $i--) {
		next unless $_[$i];
		my $d = log(abs($nn * $_[$i]/$vv)) / ($nn - $i);
		$max = $d if $d > $max;
	}
	exp($max)
}

sub poly_solve {
	my (@roots);

	if (!$_[0]) {
		push @roots, 0;
		shift until $_[0]
	}

	return sort { $a <=> $b } (@roots, &quad_solve)	if @_ <= 3;

	my @d = &diff;
	my @dr = poly_solve(@d);

	DBG and do {
		print "POLY	@_\n";
		print "DRoot	@dr\n";
	};

	my ($sl, $sr) =	
		(1 & @_) ? ($_[-1] > 0)
				? ( 1, 1) : (-1,-1)
			 : ($_[-1] > 0)
				? (-1, 1) : ( 1,-1);

	my $bnd = &upper_bound;
	my @sn = [-$bnd, $sl];
	my $max_v = 0;
	for my $x (@dr) {
		my $v = $_[-1];
		for (my $i = $#_ - 1; $i >= 0; $i--) {
			$v = $v * $x + $_[$i]
		}
		push @sn, [ $x, $v <=> 0 ];
	}
	push @sn, [ $bnd, $sr ];

	if (DBG) {
		my @colors = (31, 33, 32);
		for (@sn) {
			my $color = $colors[$_->[1] + 1];
			printf("SIGN\t%14.12f \033[${color}m(%g)\033[m\n",
				$_->[0], peval(\@_, $_->[0]));
		}
		print "\n";
	}
	for (my $i = 0; $i < $#sn; $i++) {
		if (!$sn[$i][1]) {
			push @roots, $sn[$i][0];
			next;
		}

		next if ($sn[$i][1] != -$sn[$i+1][1]);
		push @roots, find_between(\@_, \@d, $sn[$i], $sn[$i+1]);
	}

	if (DBG) {
		for (@roots) {
			printf "ROOT	%16.14f\t%g\t%g\n",
				$_, peval(\@_, $_), peval(\@d, $_);
		}
		print "=" x 40, "\n\n";
	}
	return $roots[0] ? @roots : sort { $a <=> $b } @roots;
}

sub find_between {
	my ($p, $d, $l, $r, $dbl) = @_;
	my ($left, $ls, $right, $rs) = (@$l, @$r);
	my ($dx, $dy, $y, $xn);

	my $x = ($left + $right) / 2;
	my $iter = 0;

	my $idx = $#$p - 2;
	while ($right - $left > EPSILON * 10) {
		$iter ++;

		# evaluate poly and deriv
		($y, $dy) = ($p->[-1] * $x + $p->[-2], $d->[-1]);
		for (my $i = $idx; $i >= 0; $i--) {
			$y  = $y * $x + $p->[$i];
			$dy = $dy* $x + $d->[$i]
		}

		( ($y<=>0) == $ls ? $left : $right ) = $x;
		$xn = $x - ( $dx = $y/$dy );
		$dx = -$dx if $dx < 0;
		last if ($dx < 100*EPSILON || $dx < 10*EPSILON * abs($x));

		last if $iter > 10 && abs($y) < 100*EPSILON;

		if ($xn <= $left || $xn >= $right) {
			$xn = ($left + $right) / 2;
		}

		$x = $xn;
		next;
	}
	DBG && do {
		print "FOUND	$left <- $x -> $right\n";
		print "	initial range: $l->[0] $r->[0]\n"
	};
	$x;
}

sub pmul {
	my $p = [ @{+shift} ];
	while (my $q = shift) {
		my $out;
		for (my $i = 0; $i < @$p; $i ++) {
			my $x = $p->[$i];
			for (my $j = 0; $j < @$q; $j ++) {
				$out->[$i+$j] += $x * $q->[$j]
			}
		}
		$p = $out;
	}
	$p
}

sub quad_solve {
	my ($c, $b, $a) = @_;
	my $d;
	if (@_ == 2 || !$a) {
		return $b ? -$c/$b : 0;
	} else {
		# should watch out for overflows...
		$b /= -2 * $a;

		$d = $b*$b - $c/$a;
		if (abs($d) < 5 * EPSILON * abs($b)) {
			# treat as machine error
			return $b, $b
		}
		return unless $d >= 0;

		$c = sqrt $d;
		return $b - $c, $b + $c
	}
}

1;

__END__

=head1 NAME

Math::Polynomial::RealRoots -- a simple-minded algorithm for finding real
roots of real-coefficient polynomials.

=head1 SYNOPSIS

    use Math::Polynomial::RealRoots 'real_roots';
    my @r = real_roots(1, 2, 3, -1);

=head1 DESCRIPTION

This modules provides one function you can import, C<real_roots>, which
takes a list of polynomial coefficients, zeroth order first, and returns
a list of all real roots (when it works, that is).

The algorithm does not scale very well in terms of polynomial degrees.
The goal is to use it for solving 2-D cubic Bezier curve intersections,
which involves at most 9th order monovariable equations.  This module
is probably not the most suitable tool to use for anything higher than
12-ish order root finding, but for things like 9th order or below, it
is I<fast>.

=head1 BUGS

It's not been very thoroughly tested yet.  There must be bugs.

